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ABSTRACT. The time evolution of an initially coherent, sinusoidal passive-scalar distur- 
bance is considered when the wavelength q is less than the length scale of the surrounding 
isotropic turbulent flow. In 64 s direct numerical simulations a Gaussian prescription for 
the average scalar amplitude breaks down after a timescale associated with the wavenum- 
ber of the disturbance and there is a transition to a new characteristic decay. The Gaus- 
sian prescription is given by expl — lq 2 w(t)\, where a form for w(t), the Lagrangian mean 
square displacement of a single fluid particle, is proposed. After the transition the de- 
cay is given by txp\ — t/rj, where r is the new characteristic timescale. If q > kt, then 
1/r = l/Tp + ] /rjfc, where fc* is the Kolmogorov wavenumber, Tp is the diffusive timescale 
and is the Kolmogorov timescale. For q < k^, there is a transition, but the new decay 
has not been characterized. 

An experiment originally proposed by de Gennes is considered in which the evolution 
of a coherent laser-induced pattern is read by a diffracting laser. The theory of this 
experiment involves the dispersion of particle pairs, but it is shown that in a certain limit 
it reduces to our single Fourier-mode problem and can be described in terms of single 
particle diffusion. The decay of a single mode after the transition in the simulation best 
describes the experiment. 


1. Introduction 

In many experimental situations an initially coherent, passive-scalar disturbance, 
6( r .0), is added to a turbulent flow. For example, a coherent disturbance can 
be approximately added as a single Fourier mode through illumination (Petit & 
Guyon 1979; Fermiger et al. 1982; Limat 1984) or by a heated screen in a wind 
tunnel (Warhaft & Lumley 1978). An example of a coherent disturbance in physical 
space is the thermal-interface experiment of Larue & Libby (1981). The subsequent 
evolution of the average scalar amplitude 6(~r * , t ), of the scalar fluctuations 6 2 (~r* ,t), 
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and scalar-velocity correlations, are of considerable interest. In the experiment of 
Larue & Libby (1981), a simple eddy diffusivity was able to account for the spread 
of the average profile. While these arguments should apply for long times over large 
length scales (when compared to the turbulent scales), it is not clear why they work 
well for disturbances within the turbulent scales. 

In this paper the problem of the decay of a single Fourier mode 

0(7\O )=9 0 cos(qx), (1) 

imposed on a decaying, isotropic, turbulent, velocity field is addressed using a mix- 
ture of simple analytic arguments and direct numerical simulation. The numerical 
method used is a 64 3 pseudospectral code with periodic boundary conditions. It 
has been shown by Herring & Kerr (1982) that for the Taylor-microscale Reynolds 
numbers R\ (see eq. 10), as low as 15, a simulation such as this will exhibit values 
for some correlations, such as the velocity-derivative skewness (see eq. 11), which 
are characteristic of fully developed turbulence. The nature of the simulation and 
its limitations are discussed in § 2. 

In § 3 we suggest that for short times the average amplitude of the decaying single 
mode is accurately given by 

-4(9,0 =< 0("r\O >= 6 0 exp\- l -q 2 w{t)\, (2) 

where t v(t) is the Lagrangian mean-square displacement of a fluid particle. For 
steady turbulence this formula reduces to w(t) ~ t, which is the standard eddy 
diffusivity equation for long times and large length scales (see Tennekes &; Lumley 
1972, p. 230). But for short times, we show that a different approximation to the 
Lagrangian mean-square displacement applies. In § 4 this estimate is compared 
with the amplitude A(q,t) of a single mode calculated by the numerical simulation. 
The comparison is good at short times, but the numerical results also suggest that 
for longer times the entire formulation breaks down and there is a transition to a 
different form of diffusion. We find that if the wavenumber, q , of the initial mode is 
less that the Kolmogorov wavenumber cutoff, kk (see eq. 12), and greater than j 
where L is the integral scale of the turbulence (see eq. 14), the transition will occur 
after a timescale associated with velocity disturbances of wavenumber q. However, if 
q is above the Kolmogorov wavenumber cutoff, the only velocity timescale available 
is the Kolmogorov timescale, Tk (see eq. 13), and the transition occurs only after 
that time, almost independent of q. The subsequent decay depends only on the 
diffusive and Kolmogorov timescales and obeys 8(t t )exp(—(t - t t )/r), where t t is 
the transition time and 1 jr — l/Vp + l/r*. Tp is the diffusive timescale q 2 D. 

In § 5 our result is compared to an experiment suggested by de Gennes (1977) 
and since performed by Guyon and others. In this experiment an initially coherent, 
passive-scalar disturbance of the form 

0(7*, 0) = 0o cos (9^)‘S'(7*) (3) 
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is created by an optical flash. The envelope function 5(T^) has a spatial extent l 
with ql >> 1. By diffraction of a laser beam from the resulting inhomogeneity in 
the index of refraction, the evolution of the intensity of this initial disturbance can 
be observed. For an idealized detector the measured quantity is 

/(g,t) =<\ j dt S('T y )exp{iq{x(t) - x(0)]}| 2 > . (4) 

As shown by de Gennes, (4) is formally equivalent to a problem in turbulent pair 
diffusion. However, if the envelope length scale, /, in (4) is large compared to 
the integral length scale of the turbulence, then the spatial integral in (4) does all 
the needed averaging, and the subsequent ensemble or time average is superfluous. 
Under these conditions, (4) reduces to the single Fourier mode problem of § 2,3 and 
comparisons between our results and the experiments are possible. 

The experimental method has been developed by Petit and Guyon (1979), and 
extended by Fermiger et al.(l982) and Limat (1984) for Poiseuille flow through a 
cylindrical pipe. In the experiments the Prandtl number is large and the wavenum- 
ber is above the Kolmogorov cutoff, both of which are just within the limits of our 
simulation. Significant differences between the experiment and simulation are that 
the experimental velocity field is not isotropic and that the experimental envelope 
length scale l is the order of g -1 , which implies that spatial averaging cannot be 
applied. Nevertheless, they find that the turbulent decay of the single mode is in- 
dependent of wavenumber and characterized by rjt, which is qualitatively similiar 
to our results for Pr = 4.0 and q > k k after the transition. 

2. Numerical Method 

Two numerical codes, described in Rogallo (1981) and Kerr (1985), were used 
to simulate the decay of single scalar modes. Both codes are three-dimensional 
pseudospectral codes with periodic boundary conditions. By spectral we mean that 
the fundamental variables which are stored and advanced in time are the Fourier- 
transformed velocity and scalar fields u( k ) and 0[ k ). In addition to the incom- 
pressible Navier-Stokes equations, each code can simulate any reasonable number 
of scalar equations. There are no significant differences between the two codes with 
respect to the current problem and they were used for different aspects of the prob- 
lem for historical reasons. The code described in Kerr (1985) is written in Fortran 
and used the CRAY-1S computer at NASA Ames Research Center for the R\ ss 45 
simulation. The time advancement was a third-order Runge-Kutta scheme (Wray 
1981) and the only dealiasing performed was truncation of the wavenumbers out- 
side a sphere of radius (N/2). With three scalars, 6.5 sec of cpu were used for each 
evaluation of the nonlinear terms without analysis. There are three evaluations per 
timestep and our simulation to t=0.35 required 48 timesteps or 1100 sec of cpu, 
which includes the cpu time for analysis. The eddy-turnover time (see eq. 14) for 
this simulation was r e = 0.48. 

The code described in Rogallo (1981) is written in Vectoral, and the CRAY-XMP 
computer at Ames was used for the R\ ss 15 simulation. The time advancement 
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was second-order Runge-Kutta and dealiasing is performed by using shifted grids on 
alternate evaluations and truncating outside a sphere of radius (4N/9). With three 
scalars, 4.0 sec of cpu were used per evaluation without analysis, and 800 sec of 
cpu for 80 timesteps were required for each simulation with three scalars, including 
analysis. The final time was t = 0.45 and the eddy-turnover time was r e = 0.9. 

The governing equations are the incompressible Navier-Stokes equation for the 
velocity and the transport equation for a passive scalar. The Navier-Stokes equation 
is 

du/dt + u • Vit = — Vp / p + vV 2 u (5) 

V • u = 0 (incompressibility). 

The nonlinear term, u • Vu, can be written in several different forms which are 
computationally convenient: the conservative form, 


V • (««), 


( 6 ) 


and the rotational form, 

u) x u 4- X 7 iu 2 . (7) 

The scalar equation can also be written two ways in incompressible flow: the con- 
vective form, 

dd/d< + u-Vt? = £>V 2 rf, (8) 

and the conservative form 


dti/dt + V -iid = DV 2 d. 


(9) 


In the absence of viscosity u and diffusivity D the equations conserve two positive- 
definite quadratic invariants: the kinetic energy of turbulent fluctuations, 


E = - < UiUi >, 


and the scalar variance, 


E e =<6 2 >, 


where the domain of integration is over all space. The fundamental dimensionless 
parameter that determines the velocity statistics is the Taylor-microscale Reynolds 
number 



( 10 ) 


where U is the characteristic velocity of the turbulence 
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and A is the Taylor microscale, 

A = < ui 2 > 1/2 / < (dui/dxi) 2 > 1/2 • 

Kerr (1985) uses the rotational form of the Navier-Stokes equation (7) and the 
conservative form of the scalar equation (8). Rogallo (1981) uses the conservative 
form of the Navier-Stokes equation (6) and the convective form of the scalar equation 
(9). The minimum wavenumber of both simulations was 1. 

The kinetic energy in both simulations was initialized with a 

q 4 exp{-q 2 ) 

spectrum, with the peak at q = 4 and the phases of individual wavenumbers chosen 
randomly. The velocity field was allowed to evolve until the velocity-derivative 
skewness, 



was -0.40. This is near its steady-state value of about -0.5 (Kerr, 1983) and indicates 
that an energy cascade has begun. At this time, for the first simulation R\ = 45 
and the Kolmogorov wavenumber 

fc* = (^)"= 33 , (12) 

which are the largest permissible values if one wants to maintain good resolution 
of the velocity field on the 64 3 grid. To demonstrate how well the velocity field 
is resolved and how well the energy has cascaded to higher wavenumbers, kinetic- 
energy spectra for three times from this simulation are shown in figure 1. The 
times given are with respect to the time where S u « -0.4 . The Reynolds number, 
Kolmogorov wavenumber, and the Kolmogorov timescale 

Tk=(j)*. (13) 

when S u « —0.4, are given for each simulation in Table 1. 

In figures 2 and 3, the time evolution of the energy and eddy-turnover times are 
given for the two simulations after S u « -0.4 . These will be used in the calculation 
of f{t) by (24) and (26). An eddy-turnover time is defined as 

r e = L/U, (14) 

where 

3t r f 

L = — / k~ l E{k)dk (15) 
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is the integral length scale and U 2 = 2E/3. L is the order of the simulation size 
and the energy decays as t ~ 2 8 in both simulations. Experimentally this exponent 
is 1.3 (Warhaft & Lumley 1978). Our value is higher for two reasons. First, the 
Reynolds numbers of our simulations are very low; they are therefore almost in the 
final period of decay. But more importantly, numerical simulations have a largest 
scale which inhibits the backward flow of energy. That is, the continual pairing 
of ever-larger vortices, a two-dimensional process which is also important in three- 
dimensions, is limited by the size of our box. Therefore, there is more energy at high 
wavenumbers and larger decay exponents than in experimental flows. We would not 
expect to see realistic decay exponents unless the mesh was at least 256 3 and the 
initial spectral peak of the kinetic energy was at a wavenumber greater than 8. 

Once a steady state has been reached, we impose a single passive-scalar mode 
on top of the “turbulent” velocity field. This is similiar to what was done ex- 
perimentally by Petit & Guyon (1979) when they illuminated their turbulent pipe 
with a pattern. The addition of the passive-scalar imposes additional resolution 
requirements. The relevant cutoff for Prandtl numbers > 1 is the Batchelor cutoff 
(Batchelor 1959), 

*B = (^) 4 . ( 16 ) 

In addition to the Kolmogorov wavenumber cutoff being less than the maximum 
wavenumber of the simulation, the Batchelor cutoff must also be less than this limit. 
For Pr = 0.7, this imposes no new limits on the Reynolds number. But to simulate 
higher Prandtl numbers, it is an important constraint. For Pr = 4.0 and R\ = 15 
in the second simulation, kg = 38, which is slightly outside the resolution of the 
simulation. Therefore, the large Prandtl numbers used in the experiment of Petit & 
Guyon (1979) {Pr ss 100) cannot be simulated. In the experiment, the wavenumber 
of the imposed scalar mode was between the Batchelor and Kolmogorov cutoffs. 
For simulation 2, with Pr = 4.0, this means q is between 19 and 31, the largest 
wavenumber of the simulation. Any mode placed between these limits is near the 
resolution limits of the simulation, which must be taken into account in the analysis 
of simulation 2 (§ 4). 

The values of the Prandtl numbers and wavenumbers of the single modes tested 
are given in Table 1. For comparison with purely diffusive decay, the decay of 
two modes from simulation 2 with time for Pr = 0.7 and R\ — 15 are plotted 
semilogrithmically in figure 4. The wavenumbers in each case are quite large, so 
for low Reynolds number, low Prandtl number, and high wavenumber, diffusive 
decay might be expected to dominate. Despite this, both modes initially decay 
much faster than the diffusive prediction and even faster than a simple exponential, 
which might be expected if a simple eddy diffusivity were added. § 3 will propose 
a simple analytic argument to explain this anamolous decay. 

3, Single Fourier Mode 

For an initial disturbance given by (1), the time evolution can be formally written 
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as (letting 6 o = 1) 


H~r *t) = exp[iqx(t)\, (17) 

where x(t) is the x component of the Lagrangian position of a fluid particle which 
was at x at time t = 0. The amplitude at time t in the initial Fourier mode can be 
calculated from the inverse Fourier transform 

A(q,t) — (l/V) J ,t)exp\— iqx). (18) 

For homogeneous turbulence the volume integral in (18) averages out the fluctua- 
tions and can be replaced by a time or ensemble average. Thus (18) can be formally 
rewritten as 

A(q,t) — < exp\iq{x(t) - x(0)}] >, (19) 

where x(t) - x(0) is the Lagrangian displacement of a fluid particle between time 0 
and time t. Since the turbulent velocity field is nearly Gaussian, its time integral 
will be more nearly Gaussian and the central-limit theorem can be applied (Lumley 
1970). If [ x{t ) — x(0)j is a Gaussian random variable, (19) simplifies to 

A{q,t) = exp\-^q 2 w{t)}, (20) 

where 

w(t) =< [x(t) — x(0)j 2 > (21) 

is the Lagrangian mean-square displacement. Once it is realized that (19) is the 
formal expression for the physical quantity of interest, these results are simple 
and not surprising. Relations similiar to this are common in the analysis of laser- 
doppler anemometry (Buchave et al. 1978) for time and length scales greater than 
the turbulent scales. We will show that these relations also provide a means of 
estimating w(t) for time and length scales below the turbulent scales. 

Equations (19) and (20) describe the ensemble average of (17) when all of the 
Fourier modes, except the excited one, are initially zero. Since the passive-scalar 
problem is linear, (20) can then be used along with Fourier analysis to describe the 
average behavior of any initial disturbance. The resulting prescription is to solve 
the initial-value problem for the ordinary diffusion equation, and then to make the 
substitution 

2 Dt — y u’(<) 

at the end to rescale the time behavior. 

The Lagrangian mean-square displacement w(t) can be written as 

w(t) = 2 Dt + /(<), (22) 
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where D is the coefficient of molecular diffusion, and 

f(t) = - f ds f ds' < li [lc (s),s) ■ lit [~x (s'),s') > . (23) 

3J o J o 

In (23) ^T(l?(.s), s) is the velocity of a fluid particle at time s whose position is 
7(5). Since the time domain s ^ s' makes the dominant contribution to the 
integral, (23) can be approximately rewritten in the form 

f{t) = ^J ds J ds'E{s')h{s,s '), (24) 

where E(s) is the turbulent kinetic energy per unit mass at time 5 , and the function 
h[s,s') represents the decay of correlations. For sufficiently short times f(t) is given 
approximately by 

f(t) « 2 ~E{0)t 2 . (25) 

For slightly longer times the decay of correlations can be crudely estimated if we 
know E(s) and let 

h(s,s') = exp\—[s - s')/t(s')\, (26) 

where r(s) is the eddy-turnover time for an energy-containing eddy at time s. If the 
turbulence were steady both E(s) and r(s) would be constant, and our approxima- 
tion would correspond to a simple Langevin model with / ~ t. Since the turbulence 
is decaying, the result is slightly more complicated, and f(t) is not proportional to 
t at long times. 


8 



4. Comparison to simulation 

To get the best comparison with the theory outlined above, we performed the 
highest Reynolds number simulation possible with a 64 3 grid, a decaying velocity 
field, and good resolution of the small scales. This is simulation 1. The initial scalar 
modes, q x = 4, 6, and 8, were picked so as to minimize resolution effects. When 
plotted semilogrithmically, their time evolution was qualitatively similiar to that in 
figure 4. A more illustrative manner in which to present their time decay is shown 
in figure 5. Here 

/(0 = -A'o9(^4-2£>f, (27) 

is found directly and plotted. It is seen that all three modes collapse to a single 
curve up to four times the Kolmogorov timescale, (13). In fact, the values of f(t) 
for the three modes agree to two significant figures up to two times r^. 

If .E'('S) and r(s) are taken from figures 2 and 3, combined with (24) and (26), and 
numerically integrated, f(t) can be estimated. The theoretically calculated values 
are plotted as a dashed line in figure 5. For simulation 1, f(t) can be conveniently 
fitted by the formula f(t) = 2.1 f o/3 , but there is no theoretical significance to this 
fit. Since the detailed behavior of f[i) depends on the particular energy decay in 
our simulation and a different energy decay is expected in experiments, the detailed 
behavior of f(t) will be different in a real flow. However, the procedure followed 
here should work well for experiments if the turbulence is nearly isotropic. 

At 1 k 4 Tk the calculated value for f(t) for q, = 8 breaks from the other modes 
and the theoretically predicted curve. To investigate this further we felt it was 
necessary to investigate modes where the wavenumber was the order of, or greater 
than, the Kolmogorov cutoff. This was inspired in part by the experiment of Petit 
&: Guyon (1979), whose wavenumber is always the order of, or greater than, kk • For 
simulation 1, kk ~ 33, which is the limit to the resolution of a 64 3 simulation. To 
get good resolution at the higher wavenumbers, it was necessary to reduce the Kol- 
mogorov wavenumber; this could only be done by decreasing the Reynolds number. 
Therefore, simulation 2 was performed, with R\ ^ 15. Besides calculating several 
wavenumbers, several Prandtl numbers were calculated to see what, if any, effect 
this had on the calculated value of f(t). The Reynolds number in simulation 2 is 
similiar to that used in Kerr (1981) on a 32 3 grid, but the additional resolution of 
the 64 3 grid allows higher scalar wavenumbers to be studied in the current simula- 
tion. The simulation in Kerr (1981) of the decay of single scalar modes shows decay 
similiar to simulation 1. 

Figure 6 presents f{t) for all the modes investigated in simulation 2. For short 
times, the prediction of (24) and (26) is in reasonable agreement with the directly 
computed value, although the agreement is not as good as for simulation 1. But 
after about one Kolmogorov timescale it is seen that for the higher modes there is 
a definite transition to an entirely different decay. This can also be seen for <j, = 20 
in figure 4. To illustrate this better the region of transition is enlarged in figure 7. 
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First notice that the time of transition and the subsequent decay are almost 
independent of Prandtl number for g t — 20 and g, = 12, although the transition for 
g, = 12 in figure 6 is not shown in figure 7. Recall that the Kolmogorov wavenumber 
for simulation 2 is 19. For all the modes, even those with g, > kk, the time of 
transition is never below t — Tk, which suggests that this is a lower bound for the 
transition time. It is also possible that because Pur maximum wavenumber is only 
25% greater than the Kolmogorov wavenumber the transition would occur at ever 
shorter times if higher wavenumbefs could be simulated. Since the only velocity 
timescale available at the transition is the Kolmogorov timescale, the subsequent 
decay for all the modes is compared to an exponential of the form exp( — t/rk ) for 
each mode. That is 

9 2 /(0 = f- (28) 

Tk 

Only for g, = 24 does this form seem to follow the subsequent decay closely, but 
the decay of the other modes would not be inconsistent with a decay of this form. 
The new decay is not due to the dissipation term in the scalar equation. That 
contribution has already been subtracted out in (27). 

One source for this transition could be numerical error. Very high modes near the 
resolution of the simulation are simulated and the scalar amplitude in figure 4 has 
decayed seven orders of magnitude for g, = 20. To check the size of the round-off 
error simulation 2 was repeated with g, = 16,20, and 24 for 40 timesteps without 
packing the 64-bit CRAY words into 32-bit words when doing out-of-core storage 
on disk. All of the other simulations use packing. No significant differences were 
seen without packing. 

If another source of numerical error were the cause, in particular high-wavenumber 
truncation error, strong fluctuations would be expect ed in the amplitude of the ini- 
tial mode and they would be strongest for g, = 24. All other modes with |g| « g,, 
which were initially zero, would have a similiar amplitude and would also fluctuate. 
Instead g, = 24 shows the fewest fluctuations. A simulation with g, = 27 was done 
and it did show strong fluctuations. 

Another possibility is that interference occurs when other modes feed their scalar 
variance back into the original mode. This is undoubtably occurring and is the 
probable source of the existing fluctuations about a smooth decay. But if it domi- 
nated, the sign of the original mode would be expected to change more often and the 
fluctuations would be greater. Instead, the only times that the mode changes sign 
are the few major deviations (greater than an factor of two) in the otherwise steady 
time evolution, and these last only a few timesteps before disappearing. To check 
the size of the other modes which could feed the original mode, the amplitudes of 
both the real and imaginary parts of g, = 24 are presented in figure 8. The sharp 
dips in the imaginary part represent where it changes sign. At the transition, the 
magnitude of the imaginary part is the same as that of the real part. This might 
seem large and is reflected to a small degree in the real part, but the frequency of 
the oscillations is much too large to be connected with the slow decay of the real 
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part. Even if a feedback mechanism were the dominant source of the new decay, it 
could be physically significant and might be important in explaining the effect of the 
envelope shape in the experiments. But we believe that with a large enough sample 
size, which would require a much larger simulation, the fluctuations caused by other 
modes feeding the original mode would disappear and that another explanation for 
the new decay regime exists. 

In figure 9 the mixed velocity /scalar derivative skewness 



in the direction of the initial wavevector of the scalar is plotted. Kerr (1985) has 
found that for an isotropic scalar in isotropic, homogeneous turbulence, the mixed 
skewness will be negative and about -0.5, with some variation, depending on both 
Reynolds number and Prandtl number in the range considered here. In fact, the 
mixed skewness in the directions perpendicular to the initial wavevector has this 
property in the current simulation. However, the mixed skewness in the direction of 
the initial wavevector has entirely different behavior. Initially, there is no source of 
correlation between the velocity and scalar, and the mixed skewness is zero. Next, as 
the scalar begins interacting with the velocity, the mixed skewness becomes negative, 
much as in Herring & Kerr (1982). Then it starts becoming strongly positive, 
especially for the higher wavenumber modes. Eventually the mixed skewness peaks 
and begins decreasing, presumably attaining its isotropic negative value at large 
enough times. 

The manner in which the mixed skewness changes with time appears to be corre- 
lated with the transitions in figures 4, 6, and 7. Where the mixed skewness peaks is 
roughly where the transition to the new decay is observed for each mode. This sug- 
gests a strong connection between the new decay rate and correlations containing 
both the velocity and scalar. Kerr (1985) used models from Tennekes (1968) and 
Lundgren (1982) to explain some of the correlations he calculated. In these mod- 
els the small-scale turbulent field is characterized by vortex tubes with large-scale 
stretching along the tubes and large-scale compression perpendicular to the tubes. 
Around vortices in the simulations there were, in addition, large values of the rate 
of strain that were locally induced and did not contribute to turbulence production 
or to the three-dimensional, kinetic-energy and scalar-variance cascades to large 
wavenumber. Only the large-scale components of the rate of strain contributed 
to the cascades, which is why the velocity-derivative skewness (ll) and the mixed 
velocity /scalar derivative skewness (29) were relatively small in the isotropic calcu- 
lations of Kerr (1985). The negative value of the mixed skewness was explained by 
the correlation between the strong scalar gradients and the compressive component 
of the large-scale rate of strain. 
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In this picture, to get a positive mixed skewness, as in figure 9, there must 
be scalar fluctuations correlated with the large-scale stretching, hence along the 
vortices. We propose that these scalar fluctuations eventually dominate the scalar 
field because of selective decay. In Kerr (1985), the strongest velocity gradients, and 
hence the maximum dissipation of the scalar, were perpendicular to the vortices. 
This suggests that regions of the initial scalar mode which are perpendicular to 
vortices will disappear the soonest. This is demonstrated in figure 10a. Since 
these are the fluctuations that would contribute to a negative mixed skewness, 
eventually the only fluctuations remaining are along the vortices, that is those that 
would contribute to a positive mixed skewness. These are shown in figure 10b. 
Thereafter, the dynamics of isotropic turbulence begin to dominate and the mixed 
skewness becomes more negative. 

The time it takes until the fluctuations along vortices dominate in figure 9 is 
one Kolmogorov timescale. Since this is found when the wavenumber of the initial 
mode is high and when the timescale of the decay of the fluctuations perpendicular 
to vortices is short, it implies that there must be few fluctuations along vortices, 
which implies that the vortices are highly aligned and supports Kerr (1985). After 
this the only mechanism by which the scalar can become decorrelated with the 
velocity and for there to be more decay, is for the velocity field to change and the 
vortices to move. For the smallest scales this occurs on the Kolmogorov timescale, 
which is the timescale of the decay of the scalar after the transition. 

For scalar modes with wavenumber much less that the Kolmogorov wavenum- 
ber, such as in simulation 1, the same dynamics are expected to apply, but not as 
strongly. Because the length scale of the initial scalar disturbance is much greater 
than the scale of the dissipation structures in the velocity field, the scalar fluctua- 
tions should selectively align along structures larger than the dissipation structures. 
The timescale on which this would happen would be associated with the timescale 
of these larger structures. The numerical results suggest that q 2 f{t) is always about 
10 at the transition (Table 1), or that the magnitude of the initial mode at the tran- 
sition is always nearly the same. This might provide a means of predicting when 
the transition will occur in experimental flows. 

To see how the scalar which was in the initial mode spreads into the other modes, 
the one-dimensional scalar-variance spectrum in the direction of the initial mode is 
given in figure 11 for simulation 2 and qi = 16. At later times the spectral peak 
has moved from the band containing the initial mode to lower wavenumber bands. 
This occurs well before the transition and a distinct spectral regime should also 
be noted. A k~ 5 / 3 line is drawn for comparison. The three-dimensional variance 
spectrum and the one-dimensional spectra in the other directions do not show such 
distinctive behavior. 

5. Relation to Experiment 

In applying (4) to experiment the envelope function should be evaluated at the 
time-evolved position T*(t) rather than at ~r . In addition, an actual detector 
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will see the scalar intensity in a band of wavenumbers rather than exactly at the 
initial wavenumber. Despite these oversimplifications some useful conclusions can 
be drawn. Equation (4) can be formally rewritten in the form 

I(q,t.) = J d~r* i /dr 2 5(r05(r2)<ex P I^ l( 0-x 1 + , 2 -, 2 (0)]> (30) 

In this form it is clear that the measured intensity depends on the turbulent diffusion 
of a pair of fluid particles. 

However, there is one limit in which this expression simplifies considerably. Sup- 
pose that the width l of the envelope function in any direction is large compared 
to the integral length scale L of the turbulence, and suppose further that the tur- 
bulence is homogeneous. Then the envelope function becomes irrelevant, and the 
volume integral in (4) extends over a sufficient range to average out the turbulent 
fluctuations. The ensemble average in (4) becomes superfluous and reduces to 

/(< 7 ,<) = A 2 {q,t). (31) 

where A(q,t ) is given by (19). Now the diffracted intensity is directly related to the 
single Fourier mode problem of § 3,4. The same result can be seen from (30). Under 
these conditions the main contribution to the double integral comes from trajectories 
xi(0 and x 2 (t) that are separated by a large distance compared to L. The motion 
of fluid particles 1 and 2 in (30) is then uncorrelated, and the exponential in (30) 
can be factored. The essential observation is that this factorization depends on 
the scale / of the envelope function, and not on the much smaller scale q~ l of the 
periodic pattern. 

If / becomes comparable to the integral length scale L, (30) must be looked at 
more closely. As a first approximation we continue to assume that the two-particle 
displacement 

?I2(0 = *2(0 - *2 - *i(0 + *i = <5*2(0 ~ <5*i(0 ( 32 ) 

is a Gaussian random variable. Equation (30) can then be written as 

Ifat) - J d~T x J £Tr + 2 5(7 > i)S'(T' 2 )ex P i-^g 2 u; 12 (0j (33) 

where 

u’ I2 (0 =< ki2(0] 2 >= u>i(0 + w 2 {t) - 2 < <5x 2 (0<5xj(0 > (34) 

For homogeneous turbulence, uq(<) and w 2 [t) do not depend on position, but the 
correlation term in (34) is a function of the initial separation. Equation (34) then 
becomes 


I{q,t) = exp\ — q 2 w(t)} J J d~r* id~r* 2 S [ T * \)S (T* 2 )exp\+q 2 < <5x 2 (0^*i(0 >] 
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This can be simplified further at short times. The exponent in the integrand will 
go as t 2 times the correlation in velocity between points 1 and 2 at time zero, 
but this does not allow an explicit integration. Whatever the initial separation of 
the pair of particles, they will eventually diffuse in an uncorrelated way, so that 
at long enough times the correction term in (35) will vanish. Thus we expect 
a correction for correlated pair diffusion which increases I(q,t) over the single- 
mode result at short times, but which vanishes at long times. In principle, if the 
experiments meet these conditions, the theory and simulations of § 3,4 should apply. 
In practice, experimental detectors are never ideal and corrections to our theory 
would be necessary in order to make proper comparisons. 

In Petit & Guyon (1979) and later experiments, the wavenumber of the initial 
mode is between the Kolmogorov and Batchelor cutoffs, the Prandtl number is 
high, and the time regime is between the Kolmogorov timescale and the large-scale 
eddy-turnover time. Simulation 2 is an attempt to simulate a regime like this. In 
more recent experiments, Limat (1984) finds that the decay obeys exp{-t/r) and 
that by removing the contribution of the diffusive decay a turbulent decay time r< 
can be found, where 1/r = l/Vp + l/r t and l/r# = q 2 D. That is, similiar to (28), 

1 = 9V(<). 

n 

Table 2 shows the dependence of the ratio T t jrk on Reynolds number in Limat ’s ex- 
periment along with the Kolmogorov length and timescales. There is no dependence 
on the wavelengths of the initial modes, p = 2n jq, which are given in Table 3. The 
order of magnitude of rt/r* is one, which is consistent with our second decay regime. 
If the turbulent timescale depended only on t*, then Tt/r^ should be constant. Ta- 
ble 2 show's that it is not, but its variation with Reynolds number is consistent with 
the effect of the envelope width. What should be compared to find the effect of 
the envelope is not A q/q = l~ l /q as in Limat (1984), but / _1 /fc* = pjl- That is, 
the dimension of the mixing structures must be small compared to the width of 
the envelope to minimize interference between modes, which would in turn increase 
the decay rate. Table 3 shows the experimental wavelength, the dispersion of the 
envelopes, and their width. Notice that l is constant. Therefore, since both the fun- 
damental decay and the width of the envelope are independent of wavenumber, the 
experimental decay rates should be independent of wavenumber. The only effect on 
Til 7^ from the envelope will be due to changes in ry, which decreases with Reynolds 
number. As ry decreases, the interference between modes will decrease, the decay 
will decrease, and rt/rjt will increase. Therefore, rt/rjt will increase with Reynolds 
number. From Table 2, which gives r]/l and u/Tk, we find that 

— =0.7 - 1.57. (36) 

Tk l 

For infinite Reynolds number or infinite envelope size, this predicts that r t = 0.7 r*, 
which is consistent with our calculations. Therefore we find excellent agreement 
between the experiment and our simulation. 
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It would be interesting if the experiments could show a transition similiar to the 
one in our calculations. This could be done by making measurements closer to the 
time that the initial mode was imposed, increasing the Kolmogorov timescale by 
reducing the turbulence intensity, or by decreasing the wavenumber of the initial 
mode. If the transition were observed, it could provide important new information 
on the structure of small-scale turbulence. 

6. Relation to practical problems 

Would either the initial decay of the single mode or the decay after the transition 
have an important effect on our understanding of most turbulent dispersion? Cer- 
tainly they would not have a strong effect on decay of large-scale fluctuations over 
large timescales. Even for length scales within the turbulent flow, the scalar fluctu- 
ations would decay so completely before the transition occurred that the transition 
would not be observed except in the few special cases where a single mode can be 
imposed. Besides the experiment of Petit & Guyon (1979) there might be some 
applications in laser doppler anemometry. For a sharp scalar profile in physical 
space, such as the temperature-interface problem of Larue & Libby (1981), the only 
effect might be a slightly enhanced initial turbulent diffusivity. In most practical 
applications it is difficult to separate the initial scalar profile from the large scale 
velocity field, and the regimes discussed here would not apply. 

We would like to thank M. Cloitre and E. Guyon for useful conversations concern- 
ing the relation between the calculations performed here and their laser experiment. 
We are grateful to L. Limat for sending his thesis. The work of two of the authors 
(M.N. and T.N.) was supported by the National Science Foundation under grant 
number MEA 79-02942. 
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Table 1: Velocity scales, scalar parameters (Prandtl number and wavenumber) 
and the transition parameter g, 2 /(t) for the two simulations. 


Simulation 1 

Trials 

R X = 

« k k 

= 33 

T k = 0.046 


Pr 

0. 

7 0 

.7 

0.7 



4 

6 


8 


q-f(t) 




7.8 


Simulation 2 

Trials = 

15 

k k = 

19 

T k = 0.086 

Pr 0.7 

0.7 

4.0 

4.0 

4.0 4.0 

4.0 

Qi 12 

20 

12 

16 

20 22 

24 

q^f(t) 9.6 

10.2 

9.6 

9.6 

10.2 8.4 

10.9 
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Table 2: Experimental time and length scales and appropriate ratios from Limat 
(1984). The envelope length scale l is assumed to be 275/z for all the experiments 
(see Table 3). 


Re 

2300 

3800 

6840 

n R (v) 

15 

10 

8 

n k (u) 

70 

50 

35 

T k 

3*400 

1700 

800 

T T /x k 

0.3 

0.4 

0.5 

n k /a (a = 275 u) 

0.25 

0. 18 

0.13 


Table 3. Experimental wavelengths, the dispersion of the envelopes, and the 
widths of the envelopes from Limat (1984). 


P (v) 

30 

50 

85 

100 

6 Vq 0 

6q o _1 

0.1 1 

0. 18 

0.31 

0.36 

a = P — 2 (u) 

q o 

272 

277 

274 

277 
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Figure 1: Kinetic-energy spectra for three times from simulation 1, a 64 3 simula- 
tion of decaying, isotropic turbulence whose initial spectrum was q 4 exp(-q 2 ) with 
a peak at q = 4. t - 0 is when the velocity-derivative skewness (11) S u -0.4 and 
the scalar fluctuations are introduced. R* — 45 at t — 0. Square: t = 0.0. Circle: 
t — 0.15. Triangle: t — 0.30. 
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Figure 2: The decay of kinetic energy E with time for simulations 1 (circle) and 
2 (square). 
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TIME/TkolMOGOROV 

Figure 4: The amplitude A(q, t) of an initially coherent, passive-scalar disturbance 
with spatial dependence cos[qx) plotted as a function of time for two cases from 
simulation 2. Square: Pr = 0.7, q — 12. Circle: Pr — 0.7, q = 20. The dashed lines 
represent purely diffusive decay for the two cases. 
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TIME/T K0UV , 0G0 R0V 

Figure 5: f(t) (27) plotted as a function of the Kolmogorov normalized time for 
three cases from simulation 1 (R\ = 45). Square: Pr = 0.7, q — A. Circle: Pr — 0.7, 
q = 6. Triangle: Pr — 0.7, q = 8. The agreement of all three wavenumbers to 
t — 3 rfc indicates that the Gaussian approximation is accurate for short times. The 
dashed curve is from the Lagrangian theory of § 3 (24). The apparent power-law 
fit of the theory has no theoretical significance. Note that q — 8 diverges strongly 
from the theory after t — At^. 
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TIME/T K0LM0G 0 R 0V 

Figure 6: /(<) (27) plotted as a function of the Kolmogorov normalized time 
for seven cases from simulation 2 = 15). Square: Pr = 0.7, q = 12. Circle: 

Pr = 0.7, q = 20. Triangle: Pr = 4.0, q = 12. Plus: Pr = 4.0, q = 16. Cross: 
Pr = 4.0. q = 20. Diamond: Pr = 4.0, q = 22. Inverted triangle: Pr — 4.0, q = 24. 
The relatively poor convergence at early times is probably due to the low Reynolds 
number of the simulation. Above t = the decay diverges from the Lagrangian 
theory (24), represented by the dashed line, in each case, diverging soonest for the 
highest wavenumbers. 
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.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

time/t KOLMOGOROV 

Figure 7: From figure 6, f(t) enlarged for t > 0.8rjt. The major discontinuity for 
Pr = 0.7, q = 20 is where the value of the initial mode is temporarily negative. 
The dashed lines are exponential decays associated with the Kolmogorov timescale, 
exp(-t/T k ). 
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Figure 8: Amplitudes of the real (upper curve) and imaginary (lower curve) parts 
of the initial mode for Pr = 0.7, q = 24 as functions of the Kolmogorov normalized 
time. The major discontinuities in the imaginary part are where it changes sign. 
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pr q 

a 0.7 12 

c 0.7 20 

1 4,0 12 

+ 4.0 16 

• 4.0 20 



0 .5 1.0 1.5 2.0 2.5 3.0 3.5 

HME/TkolMOGOROV 

Figure 9: The mixed velocity /scalar derivative skewnesses S u e (29) from simula- 
tion 2 as functions of the Kolmogorov normalized time. 
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Figure 10: How scalar fluctuations can align with the voriticity. Hypothetical 
scalar fluctuations are shown by sine waves, a: Looking at a three-dimensional 
vortex from the top down. The arrow indicates a streamline. The wavenumber 
of the scalar fluctuation is perpendicular to the vorticity. b: Looking at three- 
dimensional vortex from the side. The arrow indicates vortex stretching. The 
wavenumber of the scalar fluctuation is parallel to the vorticity. 
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Figure 11: One-dimensional scalar-variance spectra Ee{k\) in the direction of the 
initial mode from simulation 2 for Pr = 4.0, q — 16. The time is with respect to 
where S u % -0.4. Square: t/r^ = 0.32. Circle: </ r fc = 0.98. Triangle: i/r*. = 1.71. 
Plus: t/rk — 2.50. The dashed line represents a k~ 5/? ‘ spectrum. 
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